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ABSTRACT 

Massive black holes are key components of the assembly and evolution of cosmic 
structures and a number of surveys are currently on-going or planned to probe the 
demographics of these objects and to gain insight into the relevant physical processes. 
Pulsar Timing Arrays (PTAs) currently provide the only means to observe gravita- 
tional radiation from massive black hole binary systems with masses ^ IO^Mq. The 
whole cosmic population produces a stochastic background that could be detectable 
with upcoming Pulsar Timing Arrays. Sources sufficiently close and/or massive gener- 
ate gravitational radiation that significantly exceeds the level of the background and 
could be individually resolved. We consider a wide range of massive black hole binary 
assembly scenarios, we investigate the distribution of the main physical parameters 
of the sources, such as masses and redshift, and explore the consequences for Pulsar 
Timing Arrays observations. Depending on the specific massive black hole population 
model, we estimate that on average at least one resolvable source produces timing 
residuals in the range '--^ 5 — 50 ns. Pulsar Timing Arrays, and in particular the future 
Square Kilometre Array (SKA), can plausibly detect these unique systems, although 
the events are likely to be rare. These observations would naturally complement on the 
high-mass end of the massive black hole distribution function future surveys carried 
out by the Laser Interferometer Space Antenna {LISA). 

Key words: black hole physics, gravitational waves - cosmology: theory - pulsars: 
general 



1 INTRODUCTION 

Massive black hole (MBH) binary systems with masses in 
the range ~ 10* — 10^° M0 are amongst the primary candi- 
date sources of gravitational waves (GWs) at ~ nHz - mHz 
frequencies (see, e.g., Haehnelt 1994; JafTe & Backer 2003; 
Wyithe & Loeb 2003, Sesana et al. 2004, Sesana et al. 2005). 
The frequency band ~ 10~^ Hz — 1 Hz will be probed by the 
Laser Interferometer Space Antenna [LISA, Bender et al. 
1998), a space-borne gravitational wave laser interferome- 
ter being developed by ESA and NASA. The observational 
window 10~® Hz — 10~® Hz is already accessible with Pul- 
sar Timing Arrays (PTAs; e.g. the Parkes radio-telescope, 
Manchester 2008). PTAs exploit the effect of GWs on the 
propagation of radio signals from a pulsar to the Earth (e.g. 
Sazhin 1978, Detweiler 1979, Bertotti et al. 1983), produc- 
ing a characteristic signature in the time of arrival (TOA) 
of radio pulses. The timing residuals of the fit of the actual 
TOA of the pulses and the TOA according to a given model. 



carry the physical information about unmodelled effects, in- 
cluding GWs (e.g. Helling & Downs 1983, Jenet et al. 2005). 
The complete Parkes PTA (Manchester 2008) , the European 
Pulsar Timing Array (Janssen et al. 2008), and NanoGrav 
are expected to improve considerably on the capabilities 
of these surveys and the planned Square Kilometer Array 
(SKA; www. skatelescope.org) will produce a major leap in 
sensitivity. 

Popular scenarios of MBH formation and evolution (e.g. 
Volonteri, Haardt & Madau 2003; Wyithe & Loeb 2003, 
Koushiappas & Zentner 2006, Malbon et al. 2007, Yoo et 
al. 2007) predict the existence of a large number of mas- 
sive black hole binaries (MBHB) emitting in the frequency 
range between ~ 10~^ Hz and 10~® Hz. PTAs can gain 
direct access to this population, and address a number of 
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unanswered questions in astrophysics (such as the assem- 
bly of galaxies and dynamical processes in galactic nuclei), 
by detecting gravitational radiation of two forms: (i) the 
stochastic GW background produced by the incoherent su- 
perposition of radiation from the whole cosmic population of 
MBHBs and (ii) GWs from individual sources that are suf- 
ficiently bright (and therefore massive and/or close) so that 
the gravitational signal stands above the root-mean-square 
(rms) value of the background. Both classes of signals are 
of great interest, and the focused effort on PTAs could lead 
to the discovery of systems difficult to detect with other 
techniques. 

The possible level of the GW background, and the con- 
sequences for observations have been explored by several 
authors (see e.g. Rajagopal & Romani 1995; Phinney 2001, 
Jaffe & Backer 2003; Jenet et al. 2005; Jenet et al. 2006; 
Sesana et al. 2008). Recently, Sesana Vecchio & Colacino 
(2008, hereinafter Paperl) studied in details the properties 
of such a signal and the astrophysical information encoded 
into it, for a comprehensive range of MBHB formation mod- 
els. As shown in Paperl, there is over a factor of 10 uncer- 
tainty in the characteristic amplitude of the MBHB gener- 
ated background in the PTA frequency window. Ifowever, 
the most optimistic estimates yield an amplitude just a fac- 
tor ~ 3 below the upper-bound placed using current data 
(Jenet et al. 2006), and near-term future observations could 
either detect such a stochastic signal or start ruling out se- 
lected MBHB population scenarios. Based on our current 
astrophysical understanding of the formation and evolution 
of MBHBs and the estimates of the sensitivity of SKA, one 
could argue that this instrument guarantees the detection of 
this signal in the frequency range 3 x 10~' Hz — 5 x 10~* Hz 
for essentially every assembly scenario that is considered at 
present. 

The background generated by the cosmic population of 
MBHBs is present across the whole observational window of 
PTAs (cf. Paperl). The Monte Carlo simulations reported in 
Paperl show clearly the presence of distinctive strong peaks 
well above the average level of the stochastic contribution 
(cf. Figure 1 and 4 in Paperl). This is to be expected, as 
individual sources can generate gravitational radiation suf- 
ficiently strong to stand above the rms value of the stochas- 
tic background. These sources are of great interest because 
they can be individually resolved and likely involve the most 
massive MBHBs in the Universe. Their observation can of- 
fer further insight into the high-mass end of the MBH(B) 
population, galaxy mergers in the low-redshift Universe and 
dynamical processes that determine the formation of MBH 
pairs and the evolution to form close binaries with orbital 
periods of the order of years. 

Some exploratory studies have been carried out about 
detecting individual signals from MBHBs in PTA data 
(Jenet et al. 2004, 2005). In this paper we study system- 
atically for a comprehensive range of assembly scenarios the 
properties, in particular the distribution of masses and red- 
shift, of the sources that give rise to detectable individual 
events; we compute the induced timing residuals and the ex- 
pected number of sources at a given timing residual level. To 
this aim, the modelling of the high-mass end of the MBHB 
population at relatively low redshift is of crucial impor- 
tance. We generate a statistically significant sample of merg- 
ing massive galaxies from the on-line Millennium database 



{http://www.g-vo.org/Millennium) and populate them with 
central MBHs according to different prescriptions (Tremaine 
et al. 2002, Mclure et al. 2006, Lauer et al. 2007, Tundo et al. 
2007). The Millennium simulation (Springel et al. 2005) cov- 
ers a comoving volume of (500//iioo)^ Mpc^ (ftioo = -ffo/100 
km s~^ Mpc"'^ is the normalized Hubble parameter), ensur- 
ing a number of massive nearby binaries adequate to con- 
struct the necessary distribution. For each model we com- 
pute the stochastic background, the expected distribution of 
bright individual sources and the value of the characteristic 
timing residual 5tg„, see Equation (|20|) . for an observation 
time T. The signal-to-noise ratio at which a source can then 
be observed scales as SNR« 5tgw/i5trms where (5trms is the 
root-mean-square level of the timing residuals noise, both 
coming from the receiver and the GW stochastic background 
contribution. In the following we summarise our main re- 
sults: 

(i) The number of detectable individual sources for dif- 
ferent thresholds of the effective induced timing residuals 
5tgw is shown in Table [T] Depending on the specific MBH 
population model, we estimate that on average at least one 
resolvable source produces timing residuals in the range 
~ 5 — 50ns. Future PTAs, and in particular SKA, can plau- 
sibly detect these unique systems; the detection is however 
by no means guaranteed, events will be rare and just above 
the detection threshold. 

(ii) As expected, the brightest signals come from very 
massive systems with > 5 x 10* Mq . Here M = 
mI^^mI^'^ /{Ml + M2Y'^ is the chirp mass of the binary 
and Ml > M2 are the two black hole masses. Most of 
the resolvable sources are located at relatively high redshift 
(0.2 < z < 1.5), and not at z ^ 1 as one would naively 
expect, giving the opportunity to probe the universe at cos- 
mological distances. 

(iii) The number of resolvable MBHBs depends on the 
actual level of the stochastic background generated by the 
whole population; here we have used the standard simpli- 
fied assumption that the background level is determined by 
having more than one source per frequency resolution bin 
of width l/T, where T is the observational time. Using this 
definition we find that at frequencies less than 10~^ Hz there 
are typically a few resolvable sources, considering T = 5 yrs, 
with residuals in the range ~ 1 uHz — 1 /xHz. As the level of 
the background decreases for increasing frequencies, fainter 
sources become visible individually. 

(iv) As a sanity check, we have compared the MBHB pop- 
ulations and stochastic background levels obtained using 
data from the Millennium simulation (adopted in this pa- 
per) with those derived by means of merger tree realisations 
based on the Extended Press & Schechter (EPS) formalism 
(considered in Paperl) and have found good agreement. This 
provides an additional validation of the results of this paper 
and Paperl. Moreover it supports that EPS merger trees, if 
handled sensibly, can offer a valuable tool for the study of 
MBH evolution even at low redshift. 

The paper is organised as follows. In Section 2 we de- 
scribe MBHB population models, in particular the range of 
scenarios considered in this paper. A short review of the tim- 
ing residuals produced by GWs generated by an individual 
binary (in circular orbit) in the data collected by PTAs is 
provided in Section 3. Section 4 contains the key results of 
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the paper: the expected timing residuals from the estimated 
population of MBHBs, including detection rates for current 
and future PTAs. We also provide a comparison between the 
stochastic background computed according to the prescrip- 
tions considered here and the results of Paperl. Summary 
and conclusions are given in Section 5. 



2 THE MASSIVE BLACK HOLE BINARY 
POPULATION 

In this section we introduce the population synthesis model 
adopted to estimate the number, and astrophysical param- 
eters of MBHBs that emit GWs in the frequency region 
probed by PTAs. The two fundamental ingredients to com- 
pute the merger rate of MBHBs are (i) the merger history 
of galaxy haloes, and (ii) the MBH population associated 
to those haloes. We discuss them in turn. In building and 
evolving the population of sources we follow exactly the 
same method as in Paperl, to which we refer the reader for 
further details, with one important difference: the galactic 
haloes merger rates are derived using the data provided by 
the Millennium simulation, and not EPS-based models. We 
will justify this choice in the next sub- Section, but note that 
the two methods yield (within the statistical error) the same 
results. This is a result that is important in itself and has 
far reaching consequences (outside the specific issues related 
to PTAs). 



2.1 From merger trees to the Millennium 
database 

In Paperl we used models based on the Extended Press & 
Schechter (EPS) formalism (Press & Schechter 1974, Lacey 
& Cole 1993, Sheth & Tormen 1999), that trace the hier- 
archical assembly of dark matter haloes through a Monte- 
Carlo approach. Although EPS-based models tend to over- 
predict the bright end of the quasar luminosity function at 
z < 1 (e.g. Marulh et al. 2006), we showed that EPS halo 
merger rates at low redshift are consistent with observations 
of close galaxy pairs (Lin et al. 2004, Bell et al. 2006, De Pro- 
pris et al. 2007). In this paper we focus on MBHBs, whose 
GWs induce timing residuals above the stochastic signal 
from the whole population, and are therefore detectable as 
individual sources. The population of low/medium-redshift 
and high-mass sources will particularly impact on the re- 
sults. At low redshift, the EPS-based merger tree outputs 
need to be handled with care. In models such as those consid- 
ered in Paperl, each realization of the Universe is obtained 
by reconstructing the merger history of about 200 dark mat- 
ter haloes, see Volonteri et al. (2003) for details. The out- 
puts of the models are a list of coalescences labelled by MBH 
masses (for a given recipe that associates a MBH mass and 
to a given dark matter halo) and redshift. These events are 
then properly weighted over the observable volume shell at 
each redshift to obtain the distribution d^N/dMdzdt (see 
Paperl), that is the coalescence rate (the number of coales- 
cences A'' per time interval dt) in the chirp mass and redshift 
interval [A4, M + dA4] and [z,z + dz], respectively. The re- 
sulting distribution is reasonably smooth over most of the 
(Ai,z) plane, but small number statistics becomes impor- 
tant at z < 0.5 and A4 > W^Mq, which is an important 



region of the parameter space when one deals with individ- 
ual sources. 

To avoid this problem, in this paper we generate distri- 
butions of coalescing MBHBs using the galaxy haloes merger 
rates derived from the on-line Millennium run database. 
The Millennium simulation (Springel et al. 2005) covers a 
volume of (500//iioo)^ Mpc^ and is the ideal tool to con- 
struct a statistically representative distribution of massive 
low/medium-redshift sources. In fact, the typical ensemble 
of events available to construct the mass function of coa- 
lescing binaries is ~ 100 times larger than in a typical EPS- 
based merger tree realization. As a first step we compile cat- 
alogues of galaxy mergers from the semi-analytical model of 
Bertone et al. (2007) applied to the Millennium run. 



2.2 Populating galaxies with massive black holes 

We need to associate to each merging galaxy in our catalogue 
a central MBH, according to some sensible prescription. The 
Bertone et al. 2007 catalogue contains many properties of 
the merging galaxies, including the bulge mass Mbuigo, and 
the bulge rest frame magnitude Mv both of the progenitors 
and of the merger remnant. This is all we need in order 
to populate a galaxy with a central MBH. The process is 
twofold. 

(i) We populate the coalescing galaxies with MBHs ac- 
cording to four different MBH-host prescriptions: 

• Mbh — Afbuige in the version given by Tundo et al. 
(2007, "Tu" models, see Table 1): 



Mbh 

M„ 



8.31 + 1.12 log 



M, 



bulge 



1011 Mo 



(1) 



with a dispersion A = 0.33 dex. 

• Mbh — Afbuigo, with a redshift dependence in the ver- 
sion given by Mclure et al. (2006, "Mc" models, see Table 
1): 



Mb 



Mbui, 



= 2.071og(l + 2)-3.09. 



(2) 



with a redshift dependent dispersion A = 0.125z + 0.25 
dex (see Figure 3 of Mclure et al., 2006). 

• Mbh — Mv as given by Lauer et al. (2007, "La" mod- 
els, see Table 1): 



M 



- 8.67 - 1.32 + 22^ 



(3) 



Mq V 2.5 

with dispersion A = 0.35 dex. 

• Mbh — (7 as given by Tremaine et al. (2002, "Tr" 
models, see Table 1): 

Mbh 



Mr; 



= 8.13 + 4.02 logcrzoo, 



(4) 



where cr2oo is the velocity dispersion in units of 200 km 
s~^, and the assumed dispersion of the relation is A = 0.3 
dex. a is obtained applying the Faber- Jackson (Faber & 
Jackson 1976) relation in the form reported by Lauer et 
al. 2007 to the values of Mv obtained by the catalogue. 

To each merging system we assign MBH masses according 
to equations ([TJ-Q so that we have the masses of the two 
MBH progenitors, Mi and AI2 . For each prescription we also 
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MbH - Mbulge 

Tundo et (2007) 


Mbh - Mbulge 
Mclure et al. (2006) 


A/bh - Mv 
Lauer et al. (2007) 


A/bh - o- 
Tremaine et al. (2002) 


Single BH accretion 


Tu-SA 


Mc-SA 


La-SA 


Tr-SA 




6.9(2.7) 


6.6(2.5) 


8.1 (3.0) 


6.2(2.5) 




1.0 


1.0 (0.8) 


1 7 (^ 9^ 
1. ( (,1.^; 


U. O \ \J.OJ 




0.2 (0.4) 


0.02 (0.1) 


0.5 (0.6) 


0.01 (0.1) 




0.04 (0.2) 


0.002 (0.04) 


0.1 (0.2) 


0.002 (0.04) 


Double BH accretion 


Tu-DA 


Mc-DA 


La-DA 


Tr-DA 




8.3(2.9) 


7.3(2.7) 


9.6(3.2) 


7.0(2.8) 




2.2(1.4) 


1.6(1.1) 


2.6(1.5) 


1.2(1.0) 




0.6(0.7) 


0.2(0.4) 


0.8(0.8) 


0.06 (0.2) 




0.2(0.4) 


0.03(0.2) 


0.3(0.5) 


0.007 (0.1) 


No accretion (before merger) 


Tu-NA 


Mc-NA 


La-NA 


Tr-NA 




6.4(2.5) 


6.0(2.4) 


6.8(2.7) 


6.0(2.5) 




1.3(1.0) 


0.5(0.6) 


1.5(1.1) 


0.5(0.6) 




0.1 (0.3) 


0.07 (0.1) 


0.1(0.3) 


0.003 (0.05) 




0.02 (0.1) 


0.001 (0.03) 


0.02 (0.1) 


-{-) 



Table 1. The table summarises the 12 models of assembly of massive black hole binary populations considered in the paper - "Tu", 
"Mc", "La" and "Tr" identify the MBH-host relation; "SA", "DA", and "NA" label the accretion mode; full details on the models are 
given in Section 2 — and the total number of individually resolvable systems N{5tgy,) for selected values of the characteristic timing 
residuals (for each model, from top to bottom <5tgw = 1,10,50 and 100 ns considering an integration time of 5 yrs). The values in the 
table are the sample mean and standard deviation within brackets computed over the 1000 Monte-Carlo realisations for each model. 



calculate the mass of the MBH remnant, Mr, using the same 
equations. In all cases (Il|-(|4]), the remnant mass is > 
Ml + M2 , consistent with the fact that MBHs are expected 
to grow predominantly by accretion. We also emphasize that 
the observed scatter is included in each relation, according 
to the observational evidence that similar bulges may host 
significantly different MBHs. 

(ii) For each MBH-host relation we consider three differ- 
ent accretion scenarios: 

• The masses of the coalescing MBHs are Mi and M2. 
That is, either no accretion occurs, and the merger rem- 
nant. Ml + M2 < Mr, sits below the predicted mass, or 
accretion is triggered after the MBHB coalescence. We la- 
bel this accretion mode as "NA" (no accretion, see Table 
1). 

Post-coalescence accretion is expected for gas-rich 
mergers, where MBH pairing and coalescence is believed 
to occur on very short timescales (Mayer et al. 2007). If 
we are to assume that during a galaxy merger the MBH 
remnant is always brought on the correct correlation with 
its host, by the combination of merging and accretion, two 
additional routes are possible. 

• Accretion is triggered before the MBHB coalescence 
and only the more massive MBH (Mi) accretes mass; in 
this case the masses of the coalescing MBHs are aMi and 
A'h, where 



/3 = 



Mr 



Ml + M2 



(6) 



Mr - M2 
Wl 



(5) 



We label this accretion mode as "SA" (single BH accre- 
tion, see Table 1). 

• Accretion is triggered before the MBHB coalescence 
and both MBHs are allowed to accrete the same fractional 
amount of mass; in this case the masses of the coalescing 
MBHs are jSMi and f3M2 , where 



We label this accretion mode as "DA" (double BH accre- 
tion, see Table 1). 

The "SA" and the "DA" modes are to be expected in gas- 
poor mergers, especially in non-equal-mass mergers, where 
the dynamical evolution of the binary is much slower (e.g., 
Yu 2002) than the infall timescale of the gas (e.g., Cox et al. 
2008). In a stellar environment the orbital decay of MBHBs 
is expected to be much longer than in a gaseous environment 
(e.g., Sesana et al. 2007, Dotti et al. 2006). 

The MBHB models that we consider here relay on two 
assumptions: all bulges host a MBH, and a MBHB always 
coalesces following the hosts' merger. Regarding the first 
assumption, dynamical processes such as gravitational re- 
coil and triple MBH interactions may deplete bulges from 
their central MBH. However, if one compares the mass func- 
tion of coalescing binaries obtained here to the results of 
EPS merger tree models, where both gravitational recoil and 
triple interactions are consistently taken into account, one 
finds that the two distributions (shown in the upper-left 
panel of figure[T} are in excellent agreement, within the sta- 
tistical uncertainties. This is because triple interactions are 
likely to eject the lighter MBH from the host, leaving behind 
a massive binary in the vast majority of the cases (Volon- 
teri Haardt & Madau 2003, Hoffman & Loeb 2007). Grav- 
itational recoil, on the other hand, may be effective in ex- 
pelling light MBHs from protogalaxies at high redshift, but 
has probably a negligible impact on the population of MBHs 
in the mass range of interest for PTA observations (Volonteri 
2007). The assumption that MBHBs coalesce within an Hub- 
ble time following the hosts' merger is justified by several 
recent studies of MBHB dynamics. The stellar distribution 
interacting with the binary may be efiiciently repopulated 
as a consequence of non-axisymmetric or triaxial galaxy po- 
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tentials (Merritt & Poon 2004, Berzcik et al. 2006), or by 
massive perturbers (Perets & Alexander 2008); moreover, if 
one considers post-Newtonian corrections to the binary evo- 
lution and the effects of eccentricity, one finds that the coa- 
lescence timescale is significantly reduced (Berentzen et al. 
2008). If the binary evolution is gas-driven, typical harden- 
ing timescales are expected to be shorter than 10* yr (Escala 
et al. 2005, Dotti et al. 2006). 

In summary, we build a total of 12 models (4 MBH-host 
prescriptions x 3 accretion modes). Hereinafter, we will re- 
fer to each model using the labels associated to the prescrip- 
tions employed in it. As an example, the model based on the 
Mbh — <y relation ("Tr") with single BH accretion prescrip- 
tion ("SA") will be referred to as Tr-SA, see Table 1 for a 
summary. 



2.3 Computing the coalescence distributions 

Assigning a MBH to each galaxy, we obtain a list of coa- 
lescences (labelled by MBH masses and redshift); the same 
output quantity given by the EPS-based merger trees, used 
as a starting point in Paperl. Each event in the list can be 
now properly weighted over the observable volume shell at 
each redshift to obtain the distribution cfiN/dMdzdt along 
the lines described in Paperl. 

A further technical detail has to be considered to 
smooth the coalescence distributions. The Millennium sim- 
ulation provides better statistics for constructing the mass 
function of merging objects, but the redshift sampling is 
rather poor. In fact, the Millennium database consists of 63 
snapshots of the whole simulation. The most recent ones are 
taken at z = 0, 0.02, 0.04, 0.064 and 0.089; we need, at least 
at low redshift, to spread the events over a continuum in z, 
to obtain a sensible distance distribution of the GW sources. 
To this aim, we decouple the {M,z) dependence in the dif- 
ferential mass function at z < 0.3 and re- write N / dMdzdt 
in the following form: 



d-^N 
dMdzdt 



•^{z) X dz 



d^N 
'dMdzdt 



•fiz) X TiM,t) (7) 



By doing so we are redistributing according to a given 
function ^(z) the average mass function obtained at z < 
0.3. This is justified since, as expected, the mass func- 
tion does not show any significant evolution for redshifts 
below 0.3. The dependence on z should be of the form 
'i'{z) oc Uq X dVc/dz, where na is the galaxy/MBH num- 
ber density and dVc /dz is the differential comoving volume 
shell. At such small redshifts the impact of merger activity 
on galaxy/MBH number density is negligible (of order of 
0.1/Gyr, e.g.. Wake et al. 2008; White et al. 2007; Masjedi 
et al. 2006; Bell et al. 2006); therefore, we assume no to 
be constant. On the other hand the Universe can be consid- 
ered Euclidean, so that the differential volume shell is just 
proportional to z^. We then obtain the coalescing MBHB 
distribution in the form 



d-^N 



C xz^ X T{M,t), 



dMdzdt 

where C is a normalization factor set by the condition 



(8) 



dM 



dz- 



d^N 
'dMdzdt 



= / dM 




0' 10= 10» ICnO' 10= 108 
chirp mass [M^] 



Figure 1. Mass function of coalescing MBHBs according to 
MBH-host relations reported by several authors. The thick his- 
tograms refer to the model labeled at the top of the panel and 
described in the text. Error bars are calculated assuming a Pois- 
son error in the number count of events from the coalescence cat- 
alogues that contribute to the chirp mass interval. In the top-left 
panel the thin histogram is the coalescing MBHB mass function 
predicted by the VHMhopk model studied in Paperl. 



dzCxz''xT{M,t).{9) 



The M distributions of coalescing binaries are shown in 
figure [T] for all the MBH-host prescriptions assuming accre- 
tion on K'h only (models Tr-SA, Tu-SA, Mc-SA, La-SA ). 
The top-left panel also shows the distribution obtained by a 
reference EPS merger tree model (VHMhopk; Volonteri, Sal- 
vaterra & Haardt 2006) used in Paperl. The agreement with 
the A/bh - o- prescription ("Tr") is good for M > 10* M©, 
although the statistical errors are large due to low number 
statistics. The discrepancy for M < 1O*M0 is due both to 
the resolution limit of the Millennium simulation and to the 
fact that we relate MBHs to bulges. However, the fact that 
our sample may be incomplete for M < 10* Mq has little (if 
any) impact on the results of this paper: the MBHB popu- 
lation observable with PTAs is by far dominated by sources 
with M > 10* Mq, as we have shown in Paperl, and even 
more so when we consider systems that can be individu- 
ally resolved. We will further discuss this point in Section 
4. Note that, as discussed by, e.g., Lauer et al. 2007 and 
Tundo et al. 2007, the high mass end of the population de- 
rived using the Mbh — o" relation ("Tr" models) drops very 
steeply. The drop is much faster than in all the other cases, 
that is for distributions inferred from the Mbh — Mbuigo or 
the Mbh - Mv relations ("Tu", "Mc", "La" models). This 
is because a seems to converge to a plateau in the limit of 
very massive galaxies (Lauer et al. 2007). 
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3 TIMING RESIDUALS FROM RESOLVED 
MASSIVE BLACK HOLE BINARIES 

The search for GWs using timing data exploits the effect of 
gravitational radiation on the propagation of the radio waves 
from one (or more) pulsar(s). The characteristic signature 
of GWs on the time of arrival (TOA) of radio pulses (e.g. 
Sazhin 1978, Detweiler 1979, Bertotti et al. 1983) is a hnear 
combination of the two independent GW polarisations. In 
practice, the analysis consists in computing the difference 
between the expected and actual TOA of pulses; the timing 
residuals contain information on all the effects that have not 
been included in the fit, including GWs. In this section we 
summarise the observed signal produced by GWs, following 
closely Jenet et al. (2004), to which we refer the reader for 
further details. 

The observed timing residual generated by a GW source 
described by the independent polarisation amplitudes fe+,x 



r{t) = -(l + cosM)[r-+(t)cos(2i/.) + rx(t)sin(2^),] (10) 

where t is the time at the receiver, fi is the opening angle 
between the GW source and the pulsar relative to Earth and 
tp is the GW polarisation angle. The two functions r+_x(i) 
are defined as 



Jo 
Jo 



t' — — (1 — cos /i) 

c 



(11) 

(12) 
(13) 



Note that r^'^ {t) and r!^-*^ [t] have the same functional 
form, and result from the integration of the time evolution of 
the polarisation amplitudes at different times, with a delay 
~ 3.3 X 10^ (d/lkpc)(l — cos /^) yr, where d is the distance 
of the pulsar from the Earth. For GWs propagating exactly 
along the Earth-pulsar direction, there is no effect on the 
TOAs {r{t) = for cos^i = ±1). 

From now on we will concentrate on the timing residu- 
als produced by binary systems in circular orbit. We model 
gravitational radiation at the leading quadrupole Newto- 
nian order, that is fully justified by the fact that binaries 
in the mass and frequency range considered here are far 
from the final merger; in fact, the time to coalescence is 
~615(X/lO^M0)-'^''^(//5 X 10"® Hz)-*/^ yr. The timing 
residuals (llOp can be written as (Jenet et al, 2004): 



r{t)^ 
where 



(14) 



(t) = a{t) [a+{l + cos^ t) cos$(t) + 2ax costsin$(t)] (15) 

In the previous expressions t is the source inclination angle, 
${t) the GW phase related to the frequency f{t) (twice the 
orbital frequency) by 



$(t) = 27r / f{t')dt 



(16) 



and 

a[,/(t)] 
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100 Mpc 
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-1/3 



(17) 



is an overall scale factor that sets the size of the residuals; D 
is the luminosity distance to the GW source. The expression 
for r^P\t) can be simply obtained from Equations p5|l , ([16} 
and (|17|) by shifting the time t —f t — d{l — cos^)/c, see 
equation p2|) and (|13|) . The two functions a+ and ax are 
the "antenna beam patterns" that depend on the source 
location in the sky and the polarisation of the wave. 

MBHBs observable with PTAs produce a quasi- 
monochromatic signal - the frequency change is ~ 3 x 
W-^iM/W"^ Mq)^^^ {f/5 X 10"* Hz)"/^ nHz yr"^ - of 
known form (though unknown parameters). The optimal 
data analysis approach to search for these signals is by means 
of the well known technique of matched-filtering. The data 
set can be represented as: 



St{t) = r{t) + Stn{t) 



(18) 



where r{t) is the contribution from GWs, and Stnit) ac- 
counts for the fluctuations due to noise; the latter contribu- 
tion is the superposition of the intrinsic noise in the measure- 
ments and the GW stochastic background from the whole 
population of MBHBs. The angle- averaged optimal signal- 
to-noise ratio at which a signal from a MBHB radiating at 
(GW) frequency ~ / can be detected using a single pulsar 
is 



SUmsif) 



(19) 



In the previous expression Stimsif) is the root-mean-square 
value of the noise level Stn at frequency /, (.) represents the 
average over the source position in the sky and orientation of 
the orbital plane, and 5tgw(/) is the characteristic amplitude 
of the timing residual over the observation time T defined 



5igw(/) = ^«(/)v/7r, 



(20) 



where the numerical pre-factor comes from the angle average 
of the amplitude of the signal: 



<^a+(l + cos^ '■)^ + 4a X cos^ 



15 



(21) 



Equation [15] is appropriate to describe observations us- 
ing a single pulsar. In reality one can take advantage of the 
several pulsars that are continuously monitored to increase 
the signal-to-noise ratio, and therefore the confidence of de- 
tection: adding coherently the residuals from several pul- 
sars - currently the Parkes PTA contains about 20 pulsars, 
and more are expected to available with future instruments 
- yields an increase in signal-to-noise ratio proportional to 
the square-root of the number of pulsars used in the observa- 
tions. We will use the characteristic amplitude of the resid- 
uals Jfgw to quantify the strength of a GW signal in PTA 
observations; St^^ can be used to compute in a straightfor- 
ward way the signal-to-noise ratio, as a function of the noise 
level and number of pulsars in the array (all of which are 
quantities that do not depend on the astrophysical model), 
and therefore asses the probability of detection of sources in 
the context of a given MBHB assembly scenario. 
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4 RESULTS 

For each of the twelve models considered in Section 2 - the 
models are the result of four MBH-host galaxy prescriptions 
and three different accretion scenarios - that encompass a 
very broad range of MBHB's assembly scenarios, we com- 
pute the number of sources that are potentially resolvable 
individually and several statistical properties of the popu- 
lation, such as the redshift, chip-mass and frequency distri- 
butions, by means of Monte-Carlo realisations of the whole 
population of MBHBs according to a given model. Before 
describing the results, we provide details about each step in 
the computation of the relevant quantities. 

The distribution given by equation ^ is straightfor- 
wardly converted into dP N / dzdMdlnfr, i.e. the comoving 
number of binaries emitting in a given logarithmic (rest- 
frame) frequency interval with chirp mass and redshift in 
the range -I- dM] and [z, z -\- dz\, along the lines de- 

scribed in Section 3 of Paperl. As in Paperl, we then assume 
that in the frequency range of interest (/ > 3 x 10~^ Hz) 
the binary evolution is driven by GW emission only. This 
is a reasonable assumption, since the coalescence timescale 
for these systems is typically ;$10® yr; any other putative 
mechanism (i.e. star ejection, gas torques) of angular mo- 
mentum removal must have an enormous efflciency to com- 
pete with radiation reaction on such short timescales. For 
each of the twelve models we estimate the GW stochastic 
background (and the corresponding rms value of the tim- 
ing residuals as a function of frequency) generated by the 
sources following the scheme described in Section 4 of Pa- 
perl. Finally we generate distribution of bright, individually 
resolvable sources by running 1000 Monte-Carlo realisations 
of the whole population of MBHBs and by selecting only 
those sources whose characteristic timing residuals, equation 
( I20p , exceed the stochastic background level. We note that 
the result of which and how many sources raise above the av- 
erage noise level depends on the duration of the observation 
T, for two reasons: (i) T affects the size of the observational 
window in frequency space, in particular the minimum fre- 
quency l/T that can be reached, and (ii) the background 
level decreases as the observation time increases (as the size 
of the frequency resolution bin A/ = l/T decreases), en- 
hancing the number of individually resolvable sources. For 
the results that are described here and summarised in table 
[T] we set T = 5 yrs; increasing the data span to 10 years, 
the background level would be slightly lower, and the ad- 
ditional resolved sources would be barely brighter than the 
background. The statistics of bright, well resolvable sources 
is basically unaffected. We can then cast the results in terms 
of the cumulative number of resolvable sources as a function 
of the timing residuals, according to: 
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(22) 



where the integral is restricted to the sources that produce 
residuals above the rms level of the stochastic background; 
if we do not consider this additional constraint, then, for 
any given Jtgw , A^((5tgw) simply returns the total number of 
sources in the Monte-Carlo realisation above that particular 
residual threshold. 

Each Monte-Carlo realisation clearly yields a different 
value for A''((5tgw) (or its distribution according to a given 
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Figure 2. Summary of the properties of the population of mas- 
sive black hole binary systems - according to model Tu-SA - 
that generate gravitational waves in the frequency window cov- 
ered by Pulsar Timing Arrays. Top-left panel: characteristic am- 
plitude of the timing residuals Stg^ (equation( [20J) as a function 
of frequency; the asterisks are the residuals generated by individ- 
ual sources and the solid line is the estimated level of the GW 
stochastic background. Top-right panel: Stgw as a function of the 
number N(Stg„) of total (dotted line) and individually resolvable 
(solid line) sources, see equation 1221 Bottom-left panel: Distri- 
bution of the number of total (dotted lines) and resolvable (solid 
lines) sources per logarithmic frequency interval dN (Stg^i) /dlog f 
as a function of the GW frequency for different values of Stgvt: 
from top to bottom 1, 10 and 100 ns, respectively. Bottom-right 
panel: distribution of the total (dotted lines) and individually 
resolvable (solid lines) number of sources per logarithmic chirp 
mass interval dN{Stgy„)/dlogM as a function of chirp mass for 
different values of itgw: from top to bottom 1, 10 and 100 ns, 
respectively. An observation time of 5 years is assumed. 



parameter); the values quoted in the next section and sum- 
marised in table [T] refer to the sample mean computed over 
the set of Monte-Carlo realisations and the sample standard 
deviation. Figure [4] quantifies the typical 1-cr error in our 
estimate of the number of sources. 



4.1 Single resolvable binaries 

The large number of Monte-Carlo realisations allows us 
to study the details of the properties of the individual 
sources in a statistical sense. We concentrate in particular 
on the physical properties of the population, such as the ex- 
pected number of sources per logarithmic frequency interval 
dN{6tg^)/d\ogf, chirp mass range dN{5tgv,)/dlogA4 and 
redshift dN{5t^) / dz, and the observable parameters, such 
as the timing residuals produced by each system and the 
overall expected number of resolvable MBHBs at a given 
level of timing residual noise. A summary of the typical 
range of information that can be extracted from the sim- 
ulations is shown in figure [2] for the specific model Tu-SA. 
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Figure 3. Left panel: The effect of the MBH-host galaxy relation, assuming that accretion always takes place for a single black hole 
before merger ("SA" models), on the number of observable systems. The plot shows the number of total (thin lines) and resolvable 
(thick lines) sources N{5tg„) as a function of 5tgw, see equation 1221 Four different MBH merger scenarios are considered: Tu-SA (solid 
line), Tr-SA (long— dashed lines), Mc-SA (short— dashed lines) and La-SA (dotted lines), see Section 2 for a description of the models. 
Right panel: The effect of the MBH accretion model on the number of individually observable systems. The plot shows the number of 
resolvable sources only, N{5tgvj) as a function of Stg^j. As reference for the MBtl-host galaxy relation, models "La" (thick lines) and "Tr" 
(thin lines) are considered. The line style is as follow: model La-SA and Tr-SA (solid lines). La- DA and Tr-DA (short— dashed lines), and 
La-NA and Tr-NA (long-dashed lines). The duration of the observation is set to T = 5 yr 



The top-left panel shows the induced residuals of each sin- 
gle source compared to the level produced by the stochas- 
tic background from the whole population; the plot clearly 
shows the importance of taking into account the additional 
"noise contribution" from the brightness of the GW sky in 
considering the detectability of resolvable systems. There 
are many sources inducing residuals above, say the 5 ns 
level, however most of them contribute to the build-up of 
the background, and are not individually resolvable. The 
expected number of bright resolvable MBHBs at frequencies 
< 10~^Hz, and at a timing level > Ins, is typically around 
ten, with fainter sources resolvable at higher frequencies. 
The top-right panel shows the mean number of individual 
sources detectable as a function of 5tg^ from 1000 Monte- 
Carlo realizations of the emitting population. In this partic- 
ular case a sensitivity of « 10 ns is required to resolve an 
individual source; for a timing precision of 100 ns there is 
a 5% chances to observe a particularly bright source. Note 
that at the 1 ns level, there are ~ 100 MBHBs contribut- 
ing to the signal, however 90% of them contribute to the 
background and only about 10 sources are individually re- 
solvable. In the two bottom panels of figure [2] we plot the 
frequency and chirp mass distributions of resolvable sources 
for selected values of the minimum detectable residual am- 
plitude (5tgw. Not surprisingly, the chirp mass of observable 
systems decreases for smaller values of the considered resid- 
ual threshold, however even for an rms level of 1 ns all the 
systems are characterised by > 10* M0. The frequency 
distribution shows instead a peak corresponding to the fre- 



quency at which the background level equals the selected 
value of 5tgw. At higher frequencies the number of sources 
drastically drops because the number of emitting binaries is 
a quite steep function of frequency, A''(/) cx f~^^^; the de- 
crease at lower frequencies is because most of the emitters 
actually contribute to the background and are not individ- 
ually resolved (as clearly shown by the dotted lines). The 
qualitative behaviour of the results obtained using different 
astrophysical models is very similar to the one described 
in figure (2] with differences that affect only the numerical 
values of the different quantities. 

The central question that we want to address in this 
paper is what is the expected number of individually resolv- 
able sources that produce an effective timing residual above 
a given value, as a function of different models of MBHB 
formation and evolution. We summarise these results in fig- 
ures |3] and ID and in Table [T] where we show the mean total 
number of individual sources that exceed a given level of 
timing residual (as defined by equation (|22|l '). as a function 
of the timing residual. The qualitative behaviour of the re- 
sults is similar for all the scenarios, but the actual numbers 
vary significantly. In fact, both MBH-host relations and ac- 
cretion prescriptions have a strong impact on the statistics 
of the bright sources and, consequently, on their detection. 
We analyse the effect of the black holes populations and of 
the accretion in turn. The left panel of figure [3] shows re- 
sults where the accretion prescription is the same ("SA"), 
but the underlying MBH-host relation changes; on the other 
hand, in the right panel, we select two MBH-host relations, 
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Figure 4. The number of individually resolvable sources for se- 
lected MBHBs assembly models. The plots show N{Stgvj) as a 
function of iStgw , see equations 1221 and 1201 for a typical model 
(Tu-SA, top panel), and, in the lower panel, the two models that 
yield the largest (La-DA: top curve) and smallest (Tr-NA: bottom 
curve) number of sources. The solid lines represent the mean value 
of N{5tgvj ) and the shaded area the 1 — tr region as computed from 
1000 Monte-Carlo realisations of each MBHB population, see also 
Table 1. 



Figure 6. Distribution of the number of resolvable sources as 
a function of the gravitational wave frequency. The plots show 
dN{Stgvj)/dlog f for different values of the characteristic ampli- 
tude of the timing residuals (from right to left 5igw = 1, 10, 50 
and 100 ns). Each panel refers to a different MBH-host relation, 
while SA mode is considered (see labels in each panel). Chirp 
mass and redshift distributions for the same models are given in 
figures [7] and \8\ 




<5 V [ns] 

Figure 7. Chirp mass distribution of the number of resolvable 
Figure 5. Top panel Cumulative mean number of resolvable sources. The plots show dN{Stgy,)/dlogM for the same charac- 

sources N{5tg^) as a function of the characteristic timing resid- teristic amplitude of the timing residuals and models as in fig- 

ual (Stgw for different mass cuts: solid line: M > 10^ Mq; dashed ures|6]and |8] 

line: 10^ < M < IO^Mq; dotted line lO'^ M0 < M < 
10^ Mq. Bottom panel: same as top panel but for different red- 
shift intervals: solid line z < 0.1; dashed line 0.1 < z < 1; dotted 
line z > I. 
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Figure 8. Redshift distribution of the resolvable sources. The 
plots show dN{Stgvj)/dlog z for the same characteristic amplitude 
of the timing residuals and models as in figures |6] and [7] 

but change the accretion history. The left panel shows that 
assuming a sensitivity threshold of 30 ns, one expects to ob- 
serve of the order of one source in the La-SA model, while 
there is only a probability « 5% for the Tr-SA model. If 
we maintain the same MBH-host relation and we consider 
different accretion scenarios for e.g. the Lauer et al. popula- 
tion, the mean number of expected sources varies by a fac- 
tor of « 5 between 0.3 (La-NA) and 1.5 (La-DA), see right 
panel. Figure [3] shows that in the most pessimistic case - Tr- 
NA model, that has the sharpest cut-off of the MBHB mass 
function bright end - a precision of « 5ns should guarantee a 
positive detection; in the optimistic La-DA case, a precision 
of « 50ns should be sufficient. In turn the timing precision 
required for positive detection is in the range 5 — 50ns, that 
is basically consistent with a factor « 10 uncertainty in the 
background level estimated in Paperl. 

The typical spread around the mean values obtained in 
the Monte-Calro realisations is shown in figure|3]for selected 
models. When N{Stgv/) 3> 1 the 1-a range is roughly the 
Poisson error around the mean (reflecting the uncorrelated 
nature of sources in each Monte Carlo realisation). When 
N{Stgvr) < 1 it can be interpreted as the probability to find 
a single source above the considered iJtgw value if the actual 
MBHB population in the Universe follows the prescription 
given by the considered model; in this case a non-detection 
is trivially consistent with the model predictions. Table 1 
provides a summary of the results for the 12 models. 

It is also interesting to investigate the mass-redshift dis- 
tribution of the detectable sources. Figure \5\ shows the ex- 
pected number of detectable individual sources (as a func- 
tion of the residuals amplitude) for different redshift and 
mass ranges. Obviously, higher timing residuals correspond 
to higher M, since the strength of the signal is proportional 
to Ai^^^, and the most likely sources to be detected fall 
in the range 10* M© < X < lO'' Mq (MBHBs with M > 



Figure 9. The characteristic amplitude of the GW stochastic 
background from the population of MBHB systems. In each panel 
the thin lines identify the estimated level of the stochastic back- 
ground assuming "SA" (solid line), "DA" (dashed line) and "NA" 
(dotted line) accretion modes. The total GW amplitude from a 
single Monte— Carlo realisation of the signal corresponding to the 
"SA" accretion mode is also shown as thick solid line. 



10 Mq produce indeed larger timing residuals, but are also 
much rarer). Interestingly, the vast majority of detectable 
sources are at redshift 0.1 ^ z ^ 1, which shows that PTAs 
could probe the medium-redshift Universe, and are unlikely 
to discover nearby sources. The reason is simply that, at 
least at small redshift, the Universe volume increases as z"^. 

A summary of the properties of individual re- 
solvable sources, dN{5tg^)/dlog f , dN {5tgw) / dlog M and 
d7V((5fgw)/dlog is given in figures |6] [7] and [S] respectively, 
for all the four MBHB population models considered here, 
with accretion limited to a single black hole prior to merger 
considering different residual thresholds St^^ = 1, 10, 50, 100 
ns. All the models show the same qualitative features, as we 
have highlighted before. The frequency distribution shown 
in figure [6] was discussed above and is the same for all the 
models. The distribution of the detectable sources as a func- 
tion of chirp mass (figure [7]) peaks at « 3 x 10* M© for all 
models assuming 5tg„ = 1 ns. Increasing 5tgw, the distribu- 
tion peak shifts towards higher M and the mean number 
of events is strongly model dependent for 6tg^ > 10 ns, cf. 
the values in Table [T] The redshift distributions (figure [8} 
consistently show a broad peak in the range 0.2 < z < 1, 
due to the volume effect previously discussed. Note that in 
the "Mc" model the peak shifts toward higher redshifts as 
5tg„ increases, because in this model similar galaxies are 
populated by more massive black holes if found at higher 
redshifts, see equation ([2]). 
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Figure 10. The characteristic amplitude of the GW stochastic 
background compared to the estimate given in Paperl. The thick 
dashed line is the stochastic background predicted by the Tr-SA 
model; the dotted lines bound the background levels computed for 
all the models investigated in this paper. For comparison, the solid 
thick line is the background predicted by the VHMhopk model 
and the shaded area the range of uncertainty of the strength of 
the signal, as reported in Paperl (see Sections 4 and 5 of Paperl 
for details). 

4.2 Stochastic background 

As a sanity check, we compare the stochastic backgrounds 
derived according to the MBH populations inferred using 
the Millennium Simulation to the predictions of EPS-based 
models reported in Paperl (the reader is referred to Paperl 
for the technical details). In figure we show Monte Carlo 
generated signals for each model; in each panel we plot the 
stochastic levels according to the three different accretion 
modes discussed in Section 2. Both the accretion prescrip- 
tion and the adopted MBH-host correlation influence the 
level of the background. If accretion occurs onto the rem- 
nant (i.e. after coalescence, "NA" models), the predicted 
characteristic amplitude of the GW background can be up 
to a factor of 3 lower with respect to models in which both 
the MBHs accrete before the final coalescence ("DA" mod- 
els); on the other hand, Mbh — cr ("Tr") models predict 
lower backgrounds compared to Mbh — Mbuige ("Tu", "Mc") 
and Mbh — My ("La") models. A comparison of these re- 
sults with those presented in Paperl is given in figure [10] At 
10~* Hz the models studied here cover a characteristic am- 
plitude range consistent with the uncertainty estimated in 
Paperl. The Tr-SA model predicts a stochastic background 
that agrees with the typical EPS-model within 30% for 
/ < 10~^Hz. The high frequency end is instead steeper; this 
effect is caused by incompleteness in the low mass end of the 
MBH population. As shown in figure[TJ the mass function of 
coalescing MBHB derived from the Millennium simulation 
is not consistent with the one obtained using the EPS for- 
malism for < 10* Mq. The weight of a single dark matter 



particle in the Millennium simulation is 8.66//1100 x 10* Mq, 
allowing the reconstruction of haloes with minimum mass of 
the order of « 5 x 10^° Mq . Assuming a barionic fraction 
of 0.1, the simulation is then incomplete for barionic struc- 
tures less massive than m 5 x 10^ Mq. We checked this by 
plotting the mass function of barionic structures and find- 
ing a sudden drop below 10^ Mq . It is then inevitable that 
in the results derived from the Millennium Simulation most 
of the MBHs with mass below a few x 10*' Mq are missing. 
Since many of these MBHs are expected to merge with more 
massive ones during cosmic history, the (spurious) lack of 
MBHs in this mass rage explains the fiattening of the mass 
function (INm /dlogM shown in the top- left panel of figure 
[1] All the backgrounds are rather similar at / > 10~^ Hz 
because all the MBH prescriptions adopted lead to similar 
MBH mass functions at Mbh < 10* Mq (this fact is inde- 
pendent of the incompleteness issue). This means that the 
slope of the background on the right of the knee has a well 
defined dependence upon the adopted MBH population: the 
more pronounced is the high mass tail of the MBH mass 
function, the steeper is the high frequency end of the GW 
background. In turn, models constructed using the Millen- 
nium simulation confirm the findings of Paperl. 



5 SUMMARY AND CONCLUSIONS 

We have investigated the ability of Pulsar Timing Arrays 
(PTAs) to resolve individual massive black hole binary sys- 
tems by detecting gravitational radiation produced during 
the in-spiral phase through its effect on the residuals of the 
time of arrival of signals from radio pulsars. We have con- 
sidered a broad range of assembly scenarios, using the data 
of the Millennium simulation to evaluate the galactic haloes 
merger rates, and a total of twelve different models that 
control the relations between the mass of the central black 
holes and the galactic haloes, and the evolution of the black 
hole masses through accretion. These models therefore cover 
qualitatively (and to large extent quantitatively) the whole 
spectrum of MBH assembly scenarios currently considered. 
Regardless of the model, we estimate that at least one re- 
solvable source is expected to produce timing residuals in 
the range ~ 5 — 50 ns, and therefore future PTAs, and in 
particular the Square Kilometre Array may be able to ob- 
serve these systems. A whistle-stop summary of the mod- 
els and results is contained in Table [T] The total number 
of visible events clearly depend on the sensitivity of PTAs 
and on the astrophysical scenario. As expected, the bright- 
est sources (for PTAs) are very massive binaries with chirp 
mass M > 5x 10* Mq. However, (initially) quite surprisingly 
most of the resolvable sources are located at relatively high 
redshift (2 > 0.2). In conjunction with the observation of 
the stochastic GW background from the whole population 
of MBHBs, the identification of individual MBHBs could 
provide new constraints on the populations of these objects 
and the relevant physical processes. 

As a by-product of the analysis, we have also estimated 
the level of the GW stochastic background produced by 
the different models, finding good agreement with the es- 
timates derived using merger tree realisations based on the 
Extended Press & Schechter formalism considered in Paperl. 
Such agreement provides a further validation of the results 
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of this paper and Paperl, it shows that we can now have 
in hand self-consistent predictions for stochastic and deter- 
ministic signals from the cosmic population of MBHBs, and 
suggests that EPS merger trees could provide a valuable ap- 
proach to the studies of MBH evolution at low-to-medium 
redshift. 

As a final word of caution, we would like to stress that 
the results of this paper clearly suffer from considerable un- 
certainties determined by the still poor quantitative infor- 
mation about several parameters that control the models. 
The spread of the predictions of the expected events is there- 
fore likely dominated by the lack of knowledge of the model 
parameters, rather than the differences between the assem- 
bly scenarios. 
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